///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_DETAIL_REDUCTIONS_H_
#define _BUOLA_MAT_DETAIL_REDUCTIONS_H_

#include <buola/mat/cexpr.h>
#include <buola/mat/detail/ccolreduction.h>
#include <buola/mat/detail/crowreduction.h>

namespace buola { namespace mat { namespace detail {

struct CReductionSum
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        OScalar<tExpr> lSum=0;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum+=pE(r,c);
        
        return lSum;
    }
};

struct CReductionProd
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        OScalar<tExpr> lSum=1;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum*=pE(r,c);
        
        return lSum;
    }
};

struct CReductionMean
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        OScalar<tExpr> lSum=0;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum+=pE(r,c);
        
        return lSum/(lRows*lCols);
    }
};

struct CReductionNorm
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        OScalar<tExpr> lSum=0;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum+=sq(pE(r,c));
        
        return sqrt(lSum);
    }
};

struct CReductionNorm2
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        OScalar<tExpr> lSum=0;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum+=sq(pE(r,c));
        
        return lSum;
    }
};

struct CReductionMin
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        ///\todo use traits for that
        OScalar<tExpr> lMin=INFINITY;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();

        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lMin=std::min(lMin,pE(r,c));
        
        return lMin;
    }
};

struct CReductionMax
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        ///\todo use traits for that
        OScalar<tExpr> lMax=-INFINITY;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        ///\todo access only once (compare speed)!
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lMax=std::max(lMax,pE(r,c));
        
        return lMax;
    }
};

struct CReductionIndMax
{
    template<typename tExpr>
    OEnableIf<QVector<tExpr>,size_t> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        ///\todo use traits for that
        OScalar<tExpr> lMax=-INFINITY;
        size_t lIndMax=-1;
        size_t lSize=pE.size();
        
        for(int i=0;i<lSize;i++)
        {
            const double &lVal=pE(i);
            if(lVal>lMax)
            {
                lMax=lVal;
                lIndMax=i;
            }
        }
        
        return lIndMax;
    }
};

struct CReductionIndMin
{
    template<typename tExpr>
    OEnableIf<QVector<tExpr>,size_t> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        ///\todo use traits for that
        OScalar<tExpr> lMin=INFINITY;
        size_t lIndMin=-1;
        size_t lSize=pE.size();
        
        for(int i=0;i<lSize;i++)
        {
            const double &lVal=pE(i);
            if(lVal<lMin)
            {
                lMin=lVal;
                lIndMin=i;
            }
        }
        
        return lIndMin;
    }
};

struct CReductionAny
{
    template<typename tExpr>
    bool operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                if(pE(r,c))
                    return true;
        
        return false;
    }
};

struct CReductionAll
{
    template<typename tExpr>
    bool operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                if(!pE(r,c))
                    return false;
        
        return true;
    }
};

struct CReductionCount
{
    template<typename tExpr>
    size_t operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        size_t lCount=0;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                if(pE(r,c))
                    lCount++;
        
        return lCount;
    }
};


struct CReductionVar
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        size_t lN=lRows*lCols;
        
        OScalar<tExpr> lMean=CReductionMean()(pE);
        
        OScalar<tExpr> lSum=0;
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum+=sq(pE(r,c)-lMean);
        
        return lSum/(lN-1);
    }
};

struct CReductionStdDev
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,detail::ACCESS_ONCE>::value,"applying reduction to non-element-accessible expression");
        
        return sqrt(CReductionVar()(pE));
    }
};

struct CReductionLogSumExp
{
    template<typename tExpr>
    OScalar<tExpr> operator()(const tExpr &pE) const
    {
        static_assert(QElementAccess<tExpr,detail::ACCESS_FAST>::value,"applying reduction to non-element-accessible expression");
        
        OScalar<tExpr> lMax=CReductionMax()(pE);

        OScalar<tExpr> lSum=0;
        size_t lRows=pE.Rows();
        size_t lCols=pE.Cols();
        
        for(int r=0;r<lRows;r++)
            for(int c=0;c<lCols;c++)
                lSum+=exp(pE(r,c)-lMax);
        return lMax+log(lSum);
    }
};

template<typename tExpr>
detail::OScalar<tExpr> logsumexp(const CExpr<tExpr> &pE)
{
    msg_info() << "called wrong logsumexp\n";
    abort();
    const auto &lE=detail::element_accessible<detail::ACCESS_FAST>(pE.Derived());
    detail::OScalar<tExpr> lMax=max(lE);
    return lMax+log(sum(lE-lMax));
}


/*namespace detail*/ } /*namespace mat*/ } /*namespace buola*/ }

#endif
